(IN-) STABILITY OF SINGULAR EQUIVARIANT SOLUTIONS TO 
THE LANDAU-LIFSHITZ-GILBERT EQUATION 

JAN BOUWE VAN DEN BERG* AND J.F. WILLIAMSt 

Abstract. In this paper we use formal asymptotic arguments to understand the stabihty proper- 
ties of equivariant solutions to the Landau-Lifshitz-Gilbert model for ferromagnets. We also analyze 
both the harmonic map heatflow and Schrodinger map flow limit cases. All asymptotic results are 
verified by detailed numerical experiments, as well as a robust topological argument. The key re- 
sult of this paper is that blowup solutions to these problems are co-dimension one and hence both 
^^ unstable and non-generic. 

Finite time blowup solutions are thus far only known to arise in the harmonic map heatfiow 
in the special case of radial symmetry. Solutions permitted to deviate from this symmetry remain 
global for all time but may, for suitable initial data, approach arbitrarily close to blowup. A careful 
asymptotic analysis of solutions near blowup shows that finite-time blowup corresponds to a saddle 

^^^ fixed point in a low dimensional dynamical system. Radial symmetry precludes motion anywhere 

^iH but on the stable manifold towards blowup. 
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^_^ The Landau-Lifshitz-Gilbert problem is not invariant under radial symmetry. Nevertheless, a 

^ , similar scenario emerges in the equivariant setting: blowup is unstable. To be more precise, blowup 

^^ is co-dimension one both within the equivariant symmetry class and in the unrestricted class of initial 

"^Vj data. The value of the parameter in the Landau-Lifshitz-Gilbert equation plays a very subdued role 

in the analysis of equivariant blowup, leading to identical blowup rates and spatial scales for all 
parameter values. One notable exception is the angle between solution in inner scale (which bubbles 
C^ off) and outer scale (which remains), which does depend on parameter values. 
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Analyzing near-blowup solutions, we find that in the inner scale these solution quickly rotate 

' ' over an angle tt. As a consequence, for the blowup solution it is natural to consider a continuation 

scenario after blowup where one immediately re-attaches a sphere (thus restoring the energy lost 

' ' in blowup), yet rotated over an angle tt. This continuation is natural since it leads to continuous 

J^ dependence on initial data. 

O 
<N 

\^ Keywords: Landau-Lifshitz-Gilbert, Harmonic map heatflow, Schrodinger map flow, 

CN| asymptotic analysis, blowup, numerical simulations, adaptive numerical methods 

^^ 1. Introduction. In this paper we are interested in the existence and stability 

^"^ of finite-time singularities of the Landau-Lifshitz-Gilbert equation for maps from the 

^^ unit disk (in the plane) to the surface of the unit sphere, m : D'^ ^ S'^: 









UTfi 

—— = am X Am — 6m. x (m.x Am.), 
dt M V y. 

m.{x,t) = m.b{x) \x\ = 1, 

, m(x,0) = mo(x). 



We will always require the damping term /3 > and take a^ + /3^ = 1 without loss 
of generality This problem preserves the length of the vector tti, i.e., |7tio(^)| = 1 for 
all X implies that \m.{x,t)\ = 1 for all positive time (for all x). 

In the case a ^ 0, P > this equation arises as a model for the exchange interaction 
between magnetic moments in a magnetic spin system on a square lattice [20l Hi] . 
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Taking a = recovers the harmonic map heatflow which is a model in nematic hquid 
crystal flow [8 . It is also of much fundamental interest in differential geometry [27 . 
Finally, the conservative case /3 = is the Schrodinger map from the disk to the 
sphere, which is a model of current study in geometry p!2l [T6l [TT] . 

Stationary solutions in all these cases are harmonic maps. This allows us to analyze 
singularity formation in a unified manner for all parameter values. Traditionally, the 
Landau-Lifshitz-Gilbert equation is posed with Neumann boundary conditions, but 
this does not affect the local structure of singularities, should they arise. We note 
that in the harmonic map literature the second term on the right-hand side of the 
differential equation ( |1.1[ ) is often rewritten using the identities 

m X {m X Am) = —Am + (Am, m)m = —Am — || Vm|pm, 

where (•,•) denotes the inner product in R^ and ||Vm|p = ^^^i ^j=i{diWj)'^^ with 
Wi the components of w. 

As discussed in much greater detail in Sections |1.1| and |1.2[ there are initial data for 



which the solution to (1.1) becomes singular in finite time. In this paper we analyze 
this blowup behaviour, in particular its stability properties under (small) perturbations 
of the initial data. Considering initial data that lead to blowup, the question is 
whether or not solutions starting from slightly different initial data also blowup. Our 
main conclusion is that blowup is an unstable co- dimension one scenario. With this 
in mind we also investigate the behavior of solutions in "near- miss" of blowup, and 
the consequences this has for the continuation of the blowup solution after its blowup 
time. 

1.1. Problem formulation. We will consider two formulations for equation 



(1.1). The first is the so-called equivariant case: using polar coordinates (r^i/j) on the 



unit disk D = D^^ these are solutions of the form 

cos{mlj)u{r^t) — sm{mlj)v{r^t) 
(1.2) m{t^r^ilj) = I sin{mlj)u{r^t) -\- cos {mlj)v{r^t) 

w{r^t) 

which have the (intertwining) symmetry property m{t,-) o R^ = R^^ o m{t,-) for 
all uj and each fixed t, where R2 is a rotation over angle uj around the origin in the 
plane R^, while R^ is a rotation over angle uj around the z-axis in R^. 

The components {u^v^w) then satisfy the pointwise constraint u^ ^ v'^ ^ w'^ = 1, as 
well the differential equations 



(1.3) 



where 



Ut = a ( vAw — ( Av ^v \ w ] -\- (3 i Au i^u + Au 

Vt = ol[ —uAw + ( Au i^u \ w\ + /3 I Av i^v + Av 



^Wt = a {uAv — vAu) + j3 {Aw + Aw) , 



^ + -^ a^d A = ii^ + ^2 ^ ^2 ^ ^ (^2 ^ y^) 
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We will take n = 1 in what follows, except in Section [6J 

Alternatively, we can parametrize the solutions on the sphere via the Euler angles: 



(1.4) 



m(t,r,V^) = 



cos['0 + (/^(r, t)] sin ^(r, t) 
sin['0 + (/^(r, t)] sin ^(r, t) 



cos 0{r,t) 
where the equations for and cp are given by 

1. sin 2(9 



(1.5) 



< 



^Ot 



P^t 



a sm t 



1 



2 \r- 
sin2l9 



•^r 



smt 



^r 



-^r 



sm 



:^r 



We note that due to the splitting ij) + ^{r^t) in (1.4) the system (1.5) has one spatial 
variable. In this equivariant case the image of one radius in the disk thus fixes the 
entire map (through rotation) and we write m{t^r) = 7Ti(t,r, 0). 

In the special case a = and /3 = 1 only, there are radially symmetric solutions of 
the form (f = constant, reducing the system to a single equation 

(1.6) 0^ = o^^^lo^ '^^^^ 
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1.2. Previous results. It is well known that not all strong solutions to the 



radially symmetric harmonic map heatflow (1.6) are global in time. Equation (1.6) is 



TT-periodic in 0. Supplemented with the (finite energy) boundary condition ^(0) = 0, 
it only has stationary solutions of the form u'^ = 2 arctan qr for g G R. Hence with 
prescribed boundary data 6>(0,t) = and 0{l^t) = 6>* > n there is no accessible 
stationary profile. However, there is an associated Lyapunov functional. 



E{t) 



L 



9r{t,r) 



sin^ 0{t, r) 



r dr 



whose only stationary points are the family u'^ . It is this paradox that leads to 
blowup: there is a finite collection of (possibly finite) times at which u{0^t) "jumps" 
from riTT to (n it l)7r, losing Air of energy |27| |6J. The structure of the local solution 
close to the jumps (in time and space) is known, which allows us to analyze the 
stability of these solutions. 

The fundamental result in this area is due to Struwe [27] who first showed that 
solutions of the harmonic map heatfiow could exhibit the type of jumps described 
above and derived what the local structure of the blowup profile is. Chen, Ding and 



Ye [TT] then used super- and sub-solution arguments, applicable only to (1.6), to show 
that finite-time blowup must occur when u{0^t) = and u{l^t) > tt. The blowup 
rate and additional structural details were determined through a careful matched 
asymptotic analysis in [29] . 

The analysis is based on the original result of Struwe who showed that any solution 
which blows up in finite time must look locally (near the blowup point) like a rescaled 
harmonic map at the so-called quasi-stationary scale. That is, there is a scale r = 
0{R{t)) on which the solution takes the form 



(1.7) 



0{t^ r) -^ 2 arctan 



Rit) 



where R{t) -^ as t -> T. 



From [29] it is known that for generic initial data one has 
(1.8) i?(t)^^_I^^ as i^T 

for some k, > and blowup time T > 0, which both depend on the initial data. This 
result is intriguing as the blowup rate is very far from the similarity rate of VT — t [3 . 



While the blowup rate (1.8) was derived in [29 for the harmonic map heatflow, i.e. 
Of = 0, /3 = 1, in this paper we demonstrate that formal asymptotics imply that 
this rate is universal for all parameter values of a and /3. Very recently, proofs of 



the blowup rate (1.8) have appeared for the harmonic map heatflow [2T as well as 
the Schrodinger map flow ^23j, i.e., the two limit cases of the Landau-Lifshitz- Gilbert 
problem. 

It is common to consider radial symmetry when analyzing the blowup dynamics of 
many react ion- diffusion equations. Typically there one can show that there must 
be blowup using radially symmetric arguments. Moreover, numerical experiments 
generically show that rescaled solutions approach radial symmetry as the blowup 
time is approached. 

For the harmonic map problem the proof of blowup solutions due to Cheng, Ding and 
Ye is completely dependent on the radial symmetry. Moreover, there are stationary 
solutions to the problem which are in the homotopy class of the initial data, but which 
are not reachable under the radial symmetry constraint. This begs the question: What 
happens when we relax the constraint of radially symmetric initial data? 

The above description is mainly restricted to the harmonic map problem {a = 0) 
which has received considerably more attention than the general Landau-Lifshitz- 
Gilbert equation. Before addressing the question of stability under non-radially sym- 
metric constraints for the harmonic map heatflow problem, we first show that blowup 
solutions are still expected for the full Landau-Lifshitz- Gilbert equation with a, /3 > 0, 
see Section |2j We note that 



E{t) =11 I Or{t, r)2 + sin^ 0{t, r) (ipr{t, r)^ + ^ j 



rdr 



is a Lyapunov functional for the equivariant problem ( |1.5[ ) as long as /3 > 0, whereas 
it is a conserved quantity for the Schrodinger map flow (/3 = 0). 

We discuss the question of stability for the full problem in a uniform manner. The 
topological argument in Section [2] suggests that blowup is co-dimension one, and this 
is indeed supported by the asymptotic analysis in Section [3] and the numerics in Sec- 
tion |4j The main quantitative and qualitative properties turn out to be independent 
of the parameter values, except for the angle between sphere that bubbles off and the 
remaining part of the solution. In Section [5] we analyze near-blowup solutions. These 
solution rotate quickly over an angle tt in the inner scale. For the blowup solution this 
implies a natural continuation scenario (leading to continuous dependence on initial 
data) after the time of blowup: the lost energy is restored immediately by re-attaching 
a sphere, rotated over an angle tt. Finally, in Section [6] we present the generalization 
to the case n > 2, followed by a succinct conclusion in Section [7J 

2. The global topological picture. We present a topological argument to 
corroborate that blowup is co-dimension one. It does not distinguish between finite 
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Fig. 2.1. One parameter families of initial data for the equivariant case; several members of 
half of each family are shown (the other half lives on the hemisphere facing away from us). Left: 
for mi) ^ N one may obtain such a family for example by stereographic projection (w.r.t. N) of 
all straight lines through the point in the plane corresponding to mij. Right: for mi, = N one can 
choose the stereographic projection of parallel lines covering the plane. 

and infinite time blowup. Since the argument relies on dissipation, it works for /3 > 0, 
but since the algebra is essentially uniform in a and /3, as we shall see in Section [3) 
we would argue that the situation for the Schrodinger map flow is the same. 

Let us first consider the equivariant case, where, as explained in Section [lT] the image 
of one radius in the disk fixes the entire map, and we write m(t, r) = m(t, r, 0). 

Let 7Ti(t,0) = N (the north pole) and r?i(t, 1) = 7715. In the notation using Euler 
angles from Section [lT] by rotational symmetry we may assume that m^ = (6^5, 0), 
Oh G [0,7r]. The only equilibrium configuration satisfying these boundary conditions 
is m = (^, ip) — (2 arctangr, 0), where q = tan(6>5/2). Note that for 6b = tt there is no 
equilibrium, hence blowup must occur for all initial data in that case [2 . 

For 6>5 G [0,7r), i.e. 7725 ^ S, we shall construct a one parameter family of initial data 
mo(r;5), and we argue that at least one of the corresponding solutions blows up. 
Since the presented argument is topological, it is robust under perturbations, hence it 
proves that blowup is (at most) co-dimension one. The matched asymptotic analysis 
in Section [3] confirms this co-dimension one nature of blowup. 

We choose one-parameter families of initial data as follows. The family of initial data 
will be parametrized by 5 G 5^, or [0, 1] with the end points identified. Let mo(r; s) 
be a continuous map from [0, 1] x S^ to 6*^, such that 

(2.1) mQ{0;s) = N and mo{l;s)=mi) for all 5. 

We may then view mo as a map from S'^ -^ S'^ by identifying {0} x S^ and {1} x S^ 



to points. Now choose any continuous family mo{r;s) satisfying (2.1) such that it 
represent a degree 1 map from S'^ to itself. One such choice is obtained by using the 
stereographic projection 



T{x,y) = 



2x 2y — l + x^ + ?/ 



1 + x^ + y^ 1 + x^ + ?/2 1 + x^ + ?/■ 



Let X5 > be such that .^^ J^^ = cos ^5, i.e. x^ = tan((7r — 6b)/2) = 1/q. Then we 
choose 



1 — r 

fno{r; s) = T I x^ -\- x^ cos(27r5) ,0:5 sin(27r5)- 



see also Figure |2.1[ For the special case that mi) = N we choose 

(2.2) fn^{r; s) = T(tan(7r(r - 1/2)), tan(7r(5 - 1/2))). 

This is just one explicit choice; any homotopy of this family of initial data that obeys 
the boundary conditions (2.1) also represents a degree 1 map on S'^. Let Xq be the 



collection of initial data obtained by taking all such homotopies. It is not hard to 
see that Xq is the space of continuous functions (with the usual supremum norm) 
satisfying the boundary conditions. Let Xi be the subset of initial data in Xq for 
which the solution to the equivariant equation ( |1.3[ ) blows up in finite or infinite time. 
The following result states that the co-dimension of Xi is at most one. In particular, 
each one parameter family of initial data that represent a degree 1 map from S'^ to 
itself has at least one member that blows up. 

Proposition 2.1. Let /3 > 0. The blowup set Xi for the equivariant equation \1.^ 
has co-dimension at most 1. 

Proof. Let mo(r, 5) be any family of initial data that, via the above identification, 
represent a degree 1 map from S'^ to itself. Let m(t, r; s) correspond to the solution 
with initial data mo(r;5). As explained above, we see from ( |2.1| ) that we may view 
mo(-; •) as a map from S'^ -^ 5^ by identifying {0} x S^ and {1} x 5^ to points. Since 
the boundary points are fixed in time, we may by the same argument view m{t,-; •) 
as a map from S"^ to itself along the entire evolution. Note that since the energy is 
a Lyapunov functional for /3 > 0, any solution tends to an equilibrium as t ^ 00. If 
none of the solutions in the family would blow up (in finite or infinite time) along the 
evolution, then all solutions converge smoothly to the unique equilibrium. In partic- 
ular, for large t the map m{t^-] ■) : S'^ -^ S'^ has its image in a small neighborhood 
of this equilibrium, hence it is contractible and thus has degree 0. Moreover, if there 
is no blowup then m(t, •; •) : S'^ ^ S'^ is continuous in t, i.e. a homotopy. This is 
clearly contradictory, and we conclude blowup must occur for at least one solution 
in the one-parameter family mo(r, s). This is a topologically robust property in the 
sense that any small perturbation of mo(r, s) also represents a degree 1 map, and the 
preceding arguments thus apply to such small perturbations of mo(r, s) as well. This 
proves that the co-dimension of Xi is at most 1. D 

One may wonder what happens when (equivariant) symmetry is lost. Although a pri- 
ori the co-dimension could be higher in that case, we will show that this is not so. 
For convenience, we only deal with boundary conditions m(t, x) = N for all x G dD^ 
which simplifies the geometric picture, but the argument can be extended to more 
general boundary conditions. 

The only equilibrium solution in this situation is m{x) = N [22 . Let rnQ{r;s) be 
the family of initial data for the equivariant case with boundary condition nib = N 



(see (2.2) and Figure 2.1). Consider now the following family of initial data for the 



-N , , -^N , , X / . , . ^ \ —AT, 



general case: 

Mq {x;s) = Mq {r^ilj;s) = I sin?/; cos?/^ \fnQ{r;s). 



We see that Mq maps Dx [0, 1] to 5'^ and Mq {dD- [0, 1]) = N, but also M^ {D'^] 0) = 
'M^{D'^]l) = TV. We may thus identify dD x [0, 1] U i:> x {0,1} to a point, and 



cos^ 


— sin^ 





sini/' 


cos^ 











1 



interpret Mq as a map from S^ to S'^. In particular, Mq represents an element in 

the homotopy group 7r3(5'^) = Z. Furthermore, upon inspection, Mq represents the 
generator of the group, since it is (a deformation of) the Hopf map (see e.g. [18 ). 
Let Xq be the collection of initial data in one parameter families obtained from all 

N 

homotopies of Mq that obey the boundary conditions 

(2.3) Mo{dD; [0, 1]) = A^, and Mq(L>2; 0) = Mo{D^; 1) = N, 

Let Xi b e th e subset of initial data in Xq for which the solution to the differential 



equation (1.1 ) blows up in finite or infinite time. As before, the co-dimension of Xi is 
at most 1, showing that dropping the equivariant symmetry does not further increase 
the instability of the blowup scenario. 

Proposition 2.2. Let (3 > and nib = N. The blowup set Xi for the general 
equation ( f^.^P has co-dimension at most 1. 



Proof. The proof is analogous to the one of Proposition 2.1, but one uses 7r3(5' ) 
instead of 7r2(5'^), i.e. the degree, to obtain the contradiction. D 

As a final remark, even though blowup is co-dimension one, this does not mean it is 
irrelevant. Clearly, by changing the initial data slightly one may avoid blowup. On 
the other hand, the arguments above indicate that blowup is caused by the topology 
of the target manifold, and one can therefore not circumvent this type of singularity 
formation by simply adding additional terms to the equation (for example a physical 
effect that works on a smaller length scale), unless additional equilibria are introduced 
which refiect the pinning of a defect. 

3. Asymptotic analysis. In this section we extend the results of [29], where 
the rate of blowup for radially symmetric solutions to the harmonic heat map problem 



( 1.6 ) was determined. We will consider both the extension to the full Landau-Lifshitz- 
Gilbert equation (i.e. a ^ 0), as well as allowing a particular class of non-radial 
perturbations. We find that blowup solutions are always unstable in the equivariant 
regime. It can be understood that the blowup solutions are separatrices between two 
distinct global behaviors. 

3.1. The inner region. We will proceed with an expansion motivated by two 
facts: (i) blowup in the harmonic map heatfiow is a quasi-static modulated stationary 
solution; (ii) the full LLG problem has the same stationary profiles as the harmonic 
map heatfiow. 

Without specifying the rescaling factor R(t) yet we introduce the rescaled variable 



m 



which transforms (1.5) to 
(3.1) 

1 sin2e/l 



^l) =13 {R^Ot - RR'ie^) + a sin 6* (i?Vt " RR'iVi) 



^« + \^^ + SS^«^« = ^ ^^'^* " ^^'^^«^ " ^ ^^'^* " ^^'^^«^ • 



inner scale 1 outer scale 1 remote scale 




j-*--r: —. \ /— 




/] 1 -^ 
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oTW)) 



0{^/T^t) 



Fig. 3.1. Regions for asymptotic analysis 



To solve this in the hmit i? ^ we pose the expansion 

i9 - 6>o + {pRR' - aR^C')Oi + 

where 



(3.2) 
(3.3) 



^0 



2arctan^, 
C(t), 



represent slow movement along the two-parameter family of equilibria = 2 arctan(r/i?) 

and (p = C. 






At the next order we have 




d^^ ^ d^ 


cos 26*0 „ .( 
^2 ^1 = ?" 


^^■^' de ' e d^ 


sin 26*0 ^6*0 d(pi 
^ sin^ 00 f^C d^ 



1. 



These equations can both be solved exactly, but we omit the algebraic details since 
for the matching we only need the asymptotic behaviour as ^ ^ oo, viz.: 



(3.6) 
(3.7) 



^i~i(lne)?'-ie', 



as ^ ^ oo, 
as ^ ^ oo. 



At this stage both R{t) and C{t) are unspecified functions. They will be determined 
through the matching of the inner (r ~ R{t)) and outer regions (r ~ \/T — t), cf. 
Figure [3lj 



3.2. The outer region. To make the mechanics of the linearization and match- 
ing as transparent as possible we shall now change variables by linearizing around the 
south pole in the formulation (|1.3[): 



(7r-6>)e^^ = u^iv, 



w 



-1. 



This recovers 

(3.8) Ut = P i Urr -\ Ur ^^^ ) -^ <^ [ "^rr -\ '^r 2 ^ 

(3.9) Vt = P Ivrr ^ -Vr ^V j - a ( ii^r + -^r ^^ 

To solve this we introduce z = u-\-iv, whence 

(3.10) zt = {f3 - ia) ( Zrr + -^r ^z ) . 



This is simply the projection of the flow from the sphere on to the tangent plane at the 
pole. Notice that in the respective limits we the recover modified linear heat equation 
(a = 0) and Schrodinger equation {j3 = 0) on the tangent plane as appropriate. 

To match the inner and outer regions we first define the similarity variables 



(3.11) 



R{t) 



, T = -ln(T-t), y = e^/V, 



To reduce confusion in what follows we shall denote 



Under this change of variables equation (3.10) becomes 

1 



(3.12) 



Zr = CZ = --Zy + (/3 - ia) ( Zyy + -Zy ^Z ) . 



A solution for this equation is 



z = ae ^^'^y 



for any a — this is just z = err in (3.10). This corresponds to the eigenfunction of the 



dominant eigenvalue of £, which governs the generic long time behaviour of solutions 



of (3.12). When we allow a to vary slowly with r, we obtain a series expansion for 



the solution of the form 



(3.13) 



z r^ e 



-r/2 



cr(r)?/ + a(r) 



ia) '2y\ny 

y 



as r ^ oo. 



We now see that the introduction of a non-zero does not affect the procedure for the 
expansion. 



Denoting a{r) = ar{r) + icJiir)^ we introduce 
(3.14) 



\r + i\i = {cFrP + CFia) + ^((7^/3 — (T^a) 

= (/3 — ia){dr + i(Ti)- 
9 



We recover 6 and (p through \z\ = it — 9, argz = y) and expand for small y: 



\^\ = e^^^/ {crry + Xr-+---] + UiV + Xi- 

r4 <j^A^ + cTiXi 



(3.15) n-0^ e--l^ J\l + Af - + ''^''Ll^ y + . . . for small y, 

\ V \J\r + A- ) 

^ GiV^^Xiy-^ ^ ... 
arg z = arctan 



(3.16) (p ~ arctan - — \ — -^ — + . . . tor small y. 

Xr A^ + A^ 4 

Here and in what follows one should be slightly careful interpreting all formulae in- 
volving the arctan because of multi-valuedness. For future reference we note that 

(3.17) argz -^ arctan — for large y, 

(3.18) arg z -^ arctan —^ = arctan t^ — arctan — for small y. 

Xr O'r P 



3.3. The matching. In order to match the inner region to the outer we first 
write the inner solution in the similarity variables: 

l9 -2 arctan Te-^/^^") 

+ e'(^RA-.«»c)(-c-"^|l„(e-"»|)+e-"'|)+... 

(3.19) ^TT-2?:^+e''/'^(m-aRc) (^+\nR- l) y + . . . , 



ip^C+ (bR^C + aRR\ e^ f ^--|- In 
(3.20) ^C-M/?C' + a|J(^+lni?-l 



6 ^y'^ ( ^ ^^'^y\ y^^ ^ 



R ) R^ 

.2 I 



y 



The matching procedure now involves setting C and R such that the expansions (3.19) 



and (3.20) agree with (3.15) and (3.16) respectively, to two orders in y: 



if : 0{y^) : C - arctan (y) , 



Oiy^): 4f/5C + «|)(^+lni?) 



1 (JiXr — (JrXi 



4 A2 + Af 
10 



To solve this we set R = e ^p{r) (with p{r) algebraic in r), and after rearranging 
terms we get 



(3.21) 



C ~ arctan 



A. 



(JrAr -\- CTiXi 



i^pc+^ip-p)) (^ -inp) - |';2VAf)''i/^ - 

Since p is defined not to change exponentially fast in r, we neglect the terms of 0{lnp). 



Using the definition of A we may simplify (3.21) to get 



(3.22) 



C ~ arctan arctan — , 



-p{p - p) ^ j3{Gr\r + Cr^Ai) + a{(Ji\r - (JrK) = CTr^^r + CTi^^i , 
—p^C ~ P{cFi\r — (TrXi) — a{arXr + CTiXi) = ai&r — (JrCTi . 



Finally, we introduce C = C + arctan ^ , so that 



(3.23) 



p-2(cr^+cr,2 



C ~ arctan 



(Jr 



-p{p -p) ^ (JrCTr + CTiCTi • 
—p O ~ (JiCr — (JrCFi • 



which is independent of a and /3. This formulation strongly suggests that the case 
Of = l,/3 = is not different from the dissipative case /3 > 0. Before solving and 
studying the system (3.23) let us recall what its solutions tell us: p{t) gives an 
algebraic correction to the blowup rate, C(r) determines the local behaviour of if 
near blowup and a describes the amplitude and orientation of the solution in self- 



similar coordinates, see (3. 17), (3. 18). In order to fully understand blowup, we need 



to solve for the blowup coordinates and determine their stability. 

The blowup solution is represented by CFiij) = C(j^(r) for some c onsta n t c G M (or 
c = oo, i.e. ar = 0), with tanC = c. In particular, equations (3.17), (3.18) show 
that there is an angle tt — arctan % between the sphere bubbling off and the solution 



remaining at/after blowup, see Figure 3.2 



By rotating the sphere we may take C = without loss of generality, i.e. ai = &i = 



and &r > (note the sign), see (3.23). The blowup dynamics is described by 
(3.24) {ar — ar)r = ar , 

and p = 2a r > 0. This equation has general solutions of the form 

1 



ar = kire^ + fe/('?") 



where /(r) 
11 



as r 



oo. 



remote scale 




Fig. 3.2. The tangent plane at the south pole can be identified with the complex plane. The 
thick curve represents z(t,-) for a time near blowup. The angle between the solution in the inner 
scale (which bubbles off) and the remote scale (which remains) is indicated. 



We can immediately set /ci = as this "instability" reflects shifts in the blowup time 
and hence is not a real instability — this is common to all blowup problems [25 . We 
note that A:2 < so that indeed a^ > as r ^ oo, and p ~ — 2A:2r^ > 0. This implies 



that dr < 0, hence argz -^ tt for large y (cf. Figure 3.2). 



At this stage we have an asymptotic description of the blowup rate and its local 
structure. Unfortunately, we do not have enough information to determine stability. 
To more carefully understand the dynamics of this system we need to linearize about 
this leading order solution to find the subsequent corrections ai^, ai^, pi and Ci in a^, 
ai^ p and C, respectively. Taking aor = /(t) ~ fe/'?", croi = ^^ Po "= 2/'(r) ~ — 2/c2/r^, 
Co = 0, we get as the system for the next order 



(3.25) 



ApiPo^PoPi 



PoPl = 4aorCrir , 

(Jor 
2poPl) = O-QrO'lr + (JirCrOr •> 



2ni 



which separates into two systems. The first one is (using po = 2aor) 

Pi = 2air , 



(3.26) 



-{Pld-Qr + O-OrPl - 2aorPl) = CTQr^lr + CTorCrir , 



which, using that aor solves ( 3.24[ ), reduces to 

{air - 0-ir)r = &ir , 

the same equation as for <Jor, and provides no additional information. The other 
system is 



(3.27) 



Ci 



r&l^Ci 



O'li 

O'liCTOr 
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O'OrO'li , 



which can be rewritten as 

r(cror^li — O-iid-Qr) = (Jii&Or — CTQr^li : 



or, again using that aor solves (3.24), 



((Jii — aii)r = &ii , 



i.e., once again equation (3.24). The asymptotic behaviour of an and C is thus given 

by (a>:i,a>:2 G M) 



T 



Ci ~ — T<Ji^ ~ — /t^i — Ai:2T e'^, 

where the exponentiahy growing terms show that blowup is unstable (the neutral 
mode corresponds to a (fixed, time-independent) rotation of the sphere). 

4. Numerical computations. To supplement the formal analysis above we 
now present some numerical experiments in the radial, equivariant and fully two- 
dimensional cases. 

4.1. Numerical methods. To reliably numerically simulate potentially singu- 
lar solutions to ( |1.1| ) one needs to use adaptivity in both time and space as well satisfy 
the constraint |m(x,t)| = 1. For the former, we use r— adaptive numerical methods 
as described in [10]. This approach is based on the moving mesh PDE approach of 
Huang and Russell [T9^ combined with scale-invariance and the Sundman transfor- 
mation in time. The expository paper [10 provides many examples of this method 
being effective for computing blowup solutions to many different problems. For the 



latter we can either use formulation (1.3) and use a projection step or regularize the 



Euler angle formulation (1.5). We have implemented both and found little difference 



in efficiency or accuracy and hence will use formulation (1.3) for all but Example 1 
as it directly follows the above asymptotic analysis. 

Full two-dimensional calculations have only been performed in the case of formulation 
( |1.1[ ) and on the unit disk. This latter fact is for numerical convenience and in no 
way affects the structure of local singularities (should they arise). Here adaptivity 
was performed using the parabolic Monge- Ampere equation as described in [9j 

4.2. Numerical results. In this Section we present a sequence of numerical 
experiments to validate the results above. For examples 1-3 we used A^ = 201 spatial 
points, the monitor function 

M = |Vm| + / \\/m\dx 

and took 

dt _ 1 
d~s~ ||M||oo 

as the rescaling between computational time s and physical time t. Example 4 was 
computed on a 61 x 61 grid using the same monitor functions. Here we took |V7ti| = 

13 



x10' 




0^ 
9 


/ / 


8 




' / 


7 






6 






X 


5 






/ 


4 






/ 


3 






/ 


2 


/ / 


1 


11 / . 


u '-'-^ 





5 


10 15 2 

r/R(s) 



Fig. 4.1. Le/t- Physical grid on which the solution was computed. Right: Computational grid 
in the region of ^ = 0. Notice there is a region of essentially constant in ^ grid trajectories in this 
region. Some trajectories are leaving this region as ^ = r/R and i? — ^ 0+ . Note, in both figures only 
every fifth grid trajectory is plotted. 





Fig. 4.2. Left: Solution on physical grid at selected times. Right: Solution on computational 
grid at same times. This clearly shows that the blowup region is very well resolved. 



^Jv^ -\- v^ -\- w'^ when using form (1.3) in radial coordinates, |V7ti| = |^r| when using 



form (1.6) and the Cartesian gradient when solving the problem in two dimensions. 
In one dimension l^c = [0^ 1] and in two dimensions Qc is the unit disk. In both cases 
the integral in the monitor function is computed in the physical variables. 

Example 1 - the Radial harmonic map 

14 




Fig. 4.3. Evolution in the rescaled spatial variable. The solutions converge to the rescaled 
arctan profile, = 2 arctan(r/2i?) as predicted (plotted under the numerical solutions). 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

t 




Fig. 4.4. Blowup of initial data (pTTp computed using\L3y. Left: Evidence of blowup. Right: 
Computational grid. Note that is very similar to Figure \41\ except that we cannot compute as far 
into the blowup. 



The first example we will consider is equation (1.6) with initial data 



(4.1) 



%{r) 



Inr. 



This case has been proven to blowup with known structure [TT] [27] and asymptotically 



calculated rate R{t). Figures 4.1 and 4.2 demonstrate the method and show how the 
adaptive scheme follows the emerging similarity structure in the underlying evolution. 
Figure |4.3| shows excellent agreement with the analytical prediction of convergence to 
the arctan profile with R{t) changing over twelve orders of magnitude. 

Example 2 - Equivariant Harmonic map 



First we reconsider the example above but using equation (1.3) with the harmonic 
map case a = 0, /3 = 1. We take the same initial as (|4.1|) and set 



/(0,t) = v{0^t) = sin(|7rr)/v2, and w{0^t) = cos(|7rr). 
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i = r/R 



Fig. 4.5. Blowup of initial data jUjjV computed using U_^. Left: 9 = arctan(Vw^ -\-v'^/w). 
Note that it again converges to the rescaled arctangent profile. The initial data is not monotone in r 
as now there is also a rotation in Lp. Right: Solutions u,v,w over time. 




0.45 0.5 



Fig. 4.6. (Left) ||V?7z||oo cls a function of j (for 7 = 0.5 the computation was stopped when 
||V?77,(-, t)||oo = le8.^ (Right) Growth and decay of ||V?77,(-, t)||oo over time for a sequence of values 
of < ^ < 0.5 (in this case the dynamics are symmetric about 7 = 1/2^. 



In Figures 4.4 and 4.5 we see the same behaviour as observed above. This is not 
surprising but a reassuring test of the numerics. 



We now consider equation ( |1.3[ ) with a = and (3 

22/ 



1 for a family of initial data 



determined via stereographic projection 

2x 



(4.2) 
where 



{u^^v^^w^) = 



1 



J/2'1 



r 



tan(— 7r/2 + ttt), and ?/ = tan(— 7r/2 + 77r), for 7 G [0,1], 



which covers the sphere as 7 varies. From the discussion of Section [2] we would expect 
blowup for a single value of 7 and decay to the stationary solution in all other cases. 
Figure 4.6 shows max(^^^) \\/m\ as a function of the parameter 7 for a sequence of 



values of 7 and initial data (4.2). 



Example 3 - Full Landau-Lifshitz- Gilbert {a > 0) We now consid er th e full 

shows 



Landau-Lifshitz- Gilbert equation with a > and (3 = vT^^c? > 0. Figure 
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4.7 



snapshots in time for a = 1/a/2 and /3 = 1/a/2 as well as max(^) \\/m\ over time 
for a sequence of values of 7 in ( |4.2[ ). There is no qualitative difference to the case 
/3 = l,a = 0. 
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Fig. 4.7. (Left) Evolution of initial data (4.2) with 7 = 0.612. . . ,a = 1/v^ and (3 = 1/^2 
(Right) Evolution of the maximum gradient for a sequence of values of j with a = (3 = 1/v^- 



Example 4 - Full Landau-Lifshitz-Gilbert (a > 0) in 2 dimensions We now 

consider the full Landau-Lifshitz-Gilbert equation with a > and /3 > 0. Figure [48 
shows snapshots in time for a = 1/2, (3 = a/3/2 with initial data (4.2 ) and 7 = 0.25 and 
a small non-radial perturbation. Figure [Xo] shows snapshots in tim e for a = 1/2, /3 = 
a/3/2 but now we have taken a larger non-radial perturbation of (4.2) and varied 7 



until we had evidence of blowup. Here ||V7ti||oo = maxj=i...3((9i7Tij)^ + {d2mj 
changes almost four orders of magnitude before the computation halts. 



N2U/2 




Fig. 4.8. Evolution of the first component mi from non-radial initial data. The (Left) Initial 
data, ||V?7z||oo = 23. (Center) ||V?77,||oo = 387 (Right) ||V?7z||oo = 12. Over time the asymmetry 
grows before the solution converges towards the radially symmetric arctan profile. 

Even though the analysis above is for radial initial data we can find solutions that 
lead to blowup with carefully tuned parameters specific to given non-radial initial 
data. This is not necessarily a true blowup solution but rather a numerical one in 
the sense that it focusses to such a degree that we cannot continue the computation. 



5. Behavior of near-blowup solutions. Instability leads to a reconfiguration 
described by a quick rotation of the sphere that had almost bubbled off. The derivation 
and asymptotics of this quick rotation are presented in Section 5.1 below. This is 
highly relevant for the problem of continuing the exceptional solution that does blowup 
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Fig. 4.9. Evolution of the first component mi from non-radial initial data. (Left) Initial data, 

||Vm||oo = 28. (Center) ||Vm||oo = 953 (Right) ||Vm||oo = 9.3e4. 



after its blowup time, as explained in Section 5.3 



5.1. The quick rotation. Here we present the asymptotics of near-blowup so- 
lutions. The inner scale (in the domain; it describes a sphere in the image, or a 
semicircle in equivariant coordinates) is given by the usual ^ = r/R{t) with 

<9 - 2 arctan ^ + {(3R'R - aR^C) [^ - ^ In ^] for large ^, 

and 

if = C{t) + {PR^C + aR'R)[]^(^^{\n(^ - 1)] for large ^. 

For the outer scale (representing a small neighborhood of the south pole S in the 
image) we introduce a fast time scale t = T ^ sH. On this time scale the dynamics 
takes place at a small spatial scale r = ef, but large compared to R{t)^ i.e., i? <C £ <C 1, 
where the solution is described by {z representing coordinates in the tangent plane at 
the south pole as in Section |3.2|) 



with solution (a = a^ + icji and 7 = 7r + ^7i) 

z = cF{t)r~^ +7(t)f + .... 
Looking at the modulus and argument of z we obtain for small r 



and 

argz = arctan \ ^ — r . 

Matching \z\ to tt — and arg2; to (/?, the matching conditions read 



f-i 



fi : -e-\/3R' - aRC') InR- "'^^ + "'^^ 



argz : f° : C ~ arctan — , 



f2 



--(/?C" + aR'R-^) InR^ ^''^^ ^f '' 
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In the remote region we have z{r) ~ qr for smah r for some <7 G C, where q ^ a ^ k2/T 
is small close to blowup, as explained in Section [3J And as before, by rotating the 
sphere we may assume that q = —go, with qo > real. By matching it follows that 

z{r) ^ —qosr for large f, hence 7^ ^ —eqo and ji ^ 0. 

Hence, by rearranging the terms we obtain 

C ~ arctan — , 



e-^R'\n{l/R) 
e-^C'R\n{l/R) 



(7^ 

—^arqoe -\- acjiq^e 
(a2+cTf)i/2 • 



Let us again remove a and /3 from the formulas by setting jij. = acFr — Pcf^ and 
fii = j3cFr + occFi. In complex notation: jir + Mi^ = (/3 + Qfz)(<Jr + cr^z). Furthermore, 
write C = C + arctan § . This leads to 

C ~ arctan ^, 

Mr 

£-2i?^ln(l/i?)- ^^^^ 



.-^i^C-ln(l/i^)^ /f' 

(Mr+Mi) / 

Looking at the matching conditions, we write fir = 2£~^R(zos C and /i^ = 2£:~^i?sin (7, 
with ^ = 0(£^), which we can neglect on this time scale. We are left with the 

dynamics of C, determined by the remaining equation 

di ~ R\n{l/R)^''' 

We see that the correct time scale is e'^ = — "^ ' ^ , which is smaller the closer we are 
to blowup (and the larger qo is). Notice that indeed e ^ R since qo = 0(1/ \nR) near 
blowup as discussed before, demonstrating self-consistent separation of spatial scales. 
The angle C thus approaches ±7r depending on the initial data, unless C = 0. The 
quick rotation due to the instability is described by ^ = sin C, and the solution, in 

the original time variable, is C{t) = ±[f H-arctan[sinh(£:~^(t — T) +co))], with cq G M. 
This shows that the blowup solution acts as a separatrix between rotations in two 
opposite directions, see Figure |5.1[ The dependence on a and f3 in this scale is only 
through the fixed rotation C = C — arctan % . 

5.2. Numerical investigation of near blowup. In the previous Section we 
saw solutions whose gradient grew dramatically and then decayed as well as those 
that show blowup. We can investigate the near blowup solutions in the context of 
the previous subsection by plotting cp = arctan('u/ii) in the region where the norm 



is large and also by plotting the dynamics in the (li, I'j-plane. Figure 5.2 shows the 
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Fig. 5.1. Left: in the (C,R) phase plane the "stable manifold'' of the blowup point acts as 
separatrix. Right: geometrically it is the boundary between a rotation over an angle tt or — tt. 





Fig. 5.2. (Left) 7 < .5. (Centre) 7 = 0.5. (Right) 7 > .5. In this sequence we plainly see the 
role of the blowup solution as a separatrix. There is blowup for 7 = 0.5 and rotation away form 
blowup in opposite directions for 7 < 0.5 than for 7 > 0.5. Here u and v have been plotted in the 
spatiotemporal regime close to blowup. 



latter for three runs with a = and /3 = 1 for three values of 7 near the critical 
value 7 = 0.5. In the two cases with 7 7^ 1/2 we see the initial motion towards the 
singularity followed by decay to a regular equilibrium whereas 7 = 1/2 leads to the 
separatrix blowup behaviour. 



5.3. Continuation after blowup. Starting from smooth initial data, solutions 
to ( |1.1| ) are unique as long as they are classical, and finite time blowup may indeed 
occur for the (radially symmetric) harmonic map heatflow pjj. At a blowup point 
(in time) the strong solution terminates (at least temporarily). On the other hand, it 
is known that weak solutions of ( |1.1| ) exist globally in time [26l [271 [II HSl ll] • It is well 
established that such weak solutions are not unique [TJ [131 [3 HH] • It is thus of interest 
to come up with criteria that select the "most appropriate" weak solution. In other 
words, how should one continue a solution after the blowup time? 



For the harmonic map heatflow in two dimensions it has been shown [26l [14] that one 
uniqueness criterion is non-increasing energy 

e(t) = i/ |Vm(t)p, 

i.e., there is exactly one weak solution that has non-increasing energy e{t) for all 
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t G [0,00). Furthermore, the energy of a solution jumps down by at least 47r at a 
singularity. 

The co-dimension one character of blowup and our analysis of near-blowup solutions 



in Section 5.1 leads us to propose a different scenario for continuation after blowup. It 
has the important advantage of continuous dependence on initial data for times after 
blowup. The scenario is identical for all parameter values a and p. Namely, consider 
an equivariant solution of ( |1.1| ) that blows up as t t T. The blowup behaviour is 
characterized by a length scale R{t) ^ as 1 1 T and 

T 

(p{r, t) ^^ and 0{r, t) ^ 2 arctan — -- for r = 0{R{t)) as 1 1 T, 

R{t) 

i.e., geometrically speaking a sphere bubbles off at t = T. Based on the analysis of 
near-blowup solutions, we propose to continue the solution for t > T by immediately 
re-attaching the sphere, rotated over an angle ir with respect to the bubbled-off sphere: 

(p{r^ t) ^^ -\-7r as 1 1 T, 

with 0(r,t) - 2 arctan ~f- for r = OiRit)), and R(t) ^ as t | T. By rotating 

R[t) 

the re-attached sphere (also referred to as a reverse bubble [28 ) over an angle tt, this 
continuation framework leads to continuous dependence on initial data, since nearby 
solution that avoid blowup undergo a rapid rotation over an angle tt, as derived in 
Section 15.11 

A solution that is continued past blowup through the re-attachment of a rotated 
sphere, does not have a monotonically decreasing energy e{t). However, the renor- 
malized energy 

^ / e(t) for t + T, 

^^ \ e(T)+47r for t = T, 

is continuous and decreases monotonically. 



In the radially symmetric harmonic map heatflow case described by (1.6) this scenario 
corresponds to 

for t < T, 
l9(0,t) = <( TT for t = T, 

27r for t > T. 

For this particular case it has been proved [30] that such a re-attachment leads to a 



unique solution for t > T. For general equivariant solutions of (1.1) such an assertion 
remains an open problem. Moreover, all conclusions in Sections [3] and [5J as they fol- 
lows from formal matched asymptotics, require mathematically rigorous justification. 

6. Other n > 2. We now summarize the calculations for n = 2, 3, . . . , following 
the same methodology as for n = 1, but now the formulas are simpler (since only the 
first term in the expansion in the outer scale is needed). As analyzed in [29 , for n > 2 
blowup is in infinite time. We shall only consider blowup with one sphere bubbling off 
(i.e. Oi G (7r,27r]) in the harmonic map flow case; the adaptation to the general case 
{a 7^ 0) is analogous to Sections p^ and p] The ^-component of the large ^ asymptotics 
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in the inner scale was already calculated in [29] : 

n \2n — 2 J 



with 



oo ^2n+l 



^n= / ,-. , .9^X9 ^^ 



(l + s2^)2 2n2sin^* 



With n > 2 equation (3.5) for cpi is now replaced by 

, (2n + l)-(2n-l)g^" ^ 

^i« + ^(iTW) ^'' ^ 

Using the boundary condition 9:^1^(0) = 0, we find 

^^« = 1 Jo (TT^^"^' 

which has asymptotic behaviour cpi^ ~ En(,'^^~^ — 2n-2 ^ as ^ ^ 00. We thus find 
that 

Since the blowup for n > 2 occurs as t ^ 00, the outer variables are just the 0{1) t 
and r (i.e. not self-similar), and the equation becomes 

Zt = {P- ia) ( Zrr + -Zr 2 ^ 

The solution is asymptotically given by (with 7 and <j complex valued) 

^ = 7(t)r^ + a(t)r-^ + . . . 

This solution needs to match into the remote region where cp = n^ 6 = tt— 2 arctan qor^^ 
with qo = tan ^^^^. Hence 7 ^ —go ^ ^• 

This leads to the matching conditions (see also Section [5?T] ) 

\z\: r-": 2R^ ^ {a^ ^ a^^^^ 



arg2; : r^ : C ~ arctan ^ , 



(Jr 



r^^: ^{PR^-^^C ^anR'R^-^^)E, "^^^^ 



As before, let us transform the equation to remove the explicit dependence on a and 
/3 by setting A^ + Xii = {f3 -\- ai){ar + (Ji). Furthermore, write C = C -\- arctan % to 
obtain 

2i?^- |A|, 

C ~ argA, 

2R^C'En - nXiqo, 
2RR'En - -A^go. 
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Hence A^ = 2R^ cos C and A^ = 2R^ sin C and the remaining system is 

En 

R' = -%-R''-^cosC, 

En 

from which we easily conclude that blowup is unstable, since the equilibria i? = and 
C = kiT are all unstable (in the C direction if k is even, and in the R direction if k is 
odd). 

7. Conclusions. In this paper we have clearly demonstrated that blowup in the 
full Landau-Lifshitz-Gilbert equation is possible but that it is not generic. Instead, we 
have identified finite-time blowup as a co-dimension one phenomenon possible only for 
specially chosen initial data. It is analogous to a saddle point along whose unstable 
manifold the fiow is much slower than on the stable one. This means that while actual 
finite time blowup occurs for initial data on a set of measure zero, there is a wide 
set of initial data for which the solution gradient does increase significantly and may 
appear to blow up in numerical simulation. 

While we agree with the results in [4l|2T] about discrete blowup in this equation, their 
computations do not indicate generic blowup in the continuous problem. Blowup in 
this problem corresponds to energy concentration at small scales and so will vanish 
on any fixed grid with limited resolution. The numerical results in [4 show changes 
in energy of little more than one order of magnitude and are resolution limited with 
h = 1/64. Instead of continuous blowup, growth in those results halt when the 
solution can no longer be resolved and the authors carefully chose to call that discrete 
blowup [4]. 

In many other problems this would not be an issue as blowup typically occurs only 
because some small scale physical effects (surface-tension, high-order diffusion, satu- 
ration, etc.) have been neglected. In those cases blowup means loss of model validity. 
However, in this problem, it is the geometry of the target manifold that leads to the 
singularity and it cannot be avoided by simply adding a regularizing term. 
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